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It is well known that conventional simulation algorithms are inefficient for the statistical descrip- 
tion of macroscopic systems exactly at the critical point due to the divergence of the corresponding 
relaxation time (critical slowing down). On the other hand the dynamics in the order parameter 
space is simplified significantly in this case due to the onset of self-similarity in the associated fluc- 
tuation patterns. As a consequence the effective action at the critical point obtains a very simple 
form. In the present work we show that this simplified action can be used in order to simulate 
efficiently the statistical properties of a macroscopic system exactly at the critical point. Using the 
proposed algorithm we generate an ensemble of configurations resembling the characteristic fractal 
geometry of the critical system related to the self-similar order parameter fluctuations. As an ex- 
ample we simulate the one-component real scalar field theory at the transition point T = Tc as a 
representative system belonging to the 3 — D Ising universality class. 
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CLc I. INTRODUCTION 

S: 

O ' The traditional simulation algorithms of statistical mechanics are proven to be inefficient for the simulation of the 

^ ,, statistical properties of macroscopic systems at a critical point. This is a consequence of the well-known mechanism 

^ ' of critical slowing down leading to a divergence of the relaxation time, independently of the algorithm used, when 

] the critical point is reached 1]. Therefore, the generation of an ensemble of configurations carrying the properties 

^ of the critical system is in practice a very difficult task as the onset of equilibrium is prerequisite. On the other 
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hand, due to the universal character of the correlations at the critical point, the dynamics of the order parameter, 
described effectively through the averaged action 0] or the constrained effective potential Q, acquire a very simple 
form 0, Q reflecting the self-similarity of the associated fluctuation pattern. The consequence of self-similarity at 
the macroscopic level is at best reflected in the fractal geometry of the domains with constant magnetization in an 
Ising ferromagnet. However, fractality in the strict mathematical sense is only deflned for an ideal critical system 
embedded in a continuous space. In real physical systems fractal geometry occurs only partially between some well 
■ defined scales. Therefore, physical fractals, contrary to the corresponding exact mathematical sets, can be defined 
I also on equidistant lattices facilitating their realization through numerical simulations. On the other hand, a realistic 
. algorithm generating critical configurations of a system at its transition point should also reproduce the (partial) 
fractal geometry of the critical clusters. There are several ways to produce sets with prescribed fractal dimension 
found in the literature [6j. The associated algorithms can either be of stochastic [7| or deterministic nature j^. 
However, although fractal geometry is strongly related to critical phenomena, there is, up to now, no algorithm 
c/3 [ capable to link directly this geometrical property with the statistical mechanics of a critical system. 
^ The present work attempts a step in this direction. The main goal is to develop an efficient algorithm for the 

^ simulation of a macroscopic system at the critical point. In fact we propose a method able to generate an ensemble 
of microstates of the considered system carrying its basic critical properties such as self-similar order parameter 
fluctuations and algebraically decaying spatial correlations |^. Using the representation of the partition function in 
terms of the effective free energy at the critical point we proceed within the saddle point approximation. Then we 
show that the associated functional measure can be expressed as a summation over piecewise constant configurations 
of the order parameter within domains (clusters) of variable size covering the entire critical system. We use these 
; ^ configurations to calculate the ensemble averaged density-density correlation function leading to a power-law form 
Ci characteristic of a fractal set. Thus we recover a relation between the fractal geometry of the critical clusters and the 
canonical ensemble of critical configurations describing the statistical properties of the system at the transition point. 

The generation of an ensemble of critical configurations is of relevance for the study of the evolution of a critical 
system under the infiuence of an external potential removing it from the initial critical state. In this case the 
determined critical ensemble is introduced as a set of initial conditions for the corresponding dynamical problem 
[To). Besides, one can also consider the time evolution of a critical system under the constraints of thermodynamical 
equilibrium so that the property of self-similarity sustains for an infinitely long time interval. This situation could be 
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actually realized in the case of biological systems for which there exist several indications that they are permanently 
on a critical state . 

The paper is organized as follows: in section II we describe the basic ingredients of the proposed algorithm for 
the construction of a random fractal measure on a lattice. In section III, subsection (A), we study in detail the one 
dimensional case. For a genuine 1 — D system with short range interactions, critical behavior is not possible |12| |. 
However, it turns out that the 1 — D version of the algorithm is optimal for the illustration of the basic steps of 
the procedure. In addition, an 1 — D effective theory describing the projection of a 3 — -D critical system, can in 
principle carry imprints of the critical state. In subsection (III B), we extend the proposed algorithm for the general 
_D— dimensional case. Finally, in section IV we summarize our concluding remarks. 



II. SADDLE POINT APPROACH FOR SINGLE CLUSTER (^-CONFIGURATIONS AT THE CRITICAL 

POINT 



There are two extreme situations where the division of the whole system to subsystems and the subsequent study of 
one of them, is a good approximation. The first is when the interactions are very weak so that statistical independence 
is valid, and the whole system can be assumed to be constituted by separated building blocks. On the opposite limit, in 
a critical system, where the correlations are very strong, we expect the emergence of self-similarity and the formation 
of fractal structures. In this case investigation of a small region offers information for the entire system. Thus, at 
the critical point, it is natural to consider the partition function Zfi of an open subdomain of a critical cluster. 
Assuming that the order parameter is an one-component real field ip, Z^i is given as: 

Zn = J 5[^]e^r[^'"l, (1) 

in terms of the effective action: 

rb,fi] = ^d^x{i(Vvp(f)f + .9/+i(f)}, (2) 

expected to be valid at the critical point T = Tc 0, 0] . In (|2J) (5 is the isothermal critical exponent 01 ■ To calculate the 
partition function Q one may use the saddle point approach developed in where it is shown that the functional 
sum in Zq is saturated by the summation over the saddle point solutions of the effective action These solutions 
encompass power-law singularities. In the sum the only significant contribution comes from those saddle points for 
which the corresponding singularities lie outside the region fl. In fact, if the distance of the location of the singularities 
from the boundary of the domain increases, the corresponding statistical weight increases too. In this case the 
saddle point solutions can be well approximated by constant functions inside Q 13]. 

Since in the present work we are mainly interested in the calculation of the density-density correlation function, it 
is necessary to extend the saddle point method to the case of the generalized effective action: 

Tci^^n]^ f d^x{h\/^{x))^+g^'+\x)+j^{x)S{x)}, (3) 
Jn ^ 

involving a source term at a; = 0. The saddle point solutions of ^ obey the evolution equation: 

V^ip{x;)~gid + l)\ipix!)f ^jSix). (4) 

Interpreting Iv^C^)! as density the associated ensemble of saddle point solutions is statistically identical to the density- 
density correlation function {\(p{x)ip{0)\), in the limit j ^ jlJl]. Note that the presence of the source at a; = does 
not affect the generality and we can use any other point as reference. 

In order to simplify the presentation of the generalized saddle point calculation, in the following we will describe 
the one-dimensional case. However, as mentioned above, critical behavior, in the absence of long range interactions, 
occurs only in higher dimensional systems. Thus, the 1 — D case should be considered as a valuable explanatory 
tool allowing also for analytical results, or as an effective description of the 1 — D projection of a higher dimensional 
critical system. In analogy with the treatment in |l3j , eq.(0J is solved as an initial value problem. The solution in 
one dimension is given implicitly in terms of the function: 
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where is a parameter depending on the assumed initial conditions and having the interpretation of the total energy 
of the system. F {a, (3;^; z) is the usual hypergeometric function |[l5|. For x > the saddle point solutions are the 
inverse function of: 

x = H{cp{x),E+)-Hiip{0),E+), (6) 

while for a; < they are the inverse function of 

x^Hiip{x),E-)-H{ipiO),E-). (7) 

Due to the source term the energy of the system varies in the two half-spaces and the corresponding energy difference 
is determined by the formula: 

E-=E+ + -^+ J ^2 +51(^(0)1^+1), (8) 

In the limit j —> 0, E- = E+ = E and these solutions simplify to those presented in [l^ . 

Similarly to the previous discussion, a typical generalized saddle point solution possesses two singularities which 
can be seen in the graphical presentation in fig^ Additionally it possesses a discontinuity in the field derivative at 




FIG. 1: A typical generalized saddle point solution, obtained analytically, for the 1 — D case. The discontinuity of the derivative 
at X = is clearly seen m the inset figure. 



X = Q which vanishes in the limit j — > 0. In analogy to the case without the source term, the summation over the 
generalized saddle points is also dominated by configurations for which the singularities lie beyond the range of the 
cluster. Therefore, the solutions I©,© which contribute significantly to the system's partition function can also be 
well approximated by a constant function within the cluster. 

The saddle point solutions in the higher dimensional case are impossible to be found analytically. However numerical 
calculations reveal a similar behavior, that is almost constant configurations with singularities forming a closed Z? — 1 
subspace around the location of the source, which can be handled as described above. Therefore, the summation over 
the saddle point solutions in the single cluster partition function, in any space dimension _D, can be well approximated 
by an ordinary integration over constant field configurations. 

The fractal geometry of the critical clusters is revealed through the power-law dependence of the density-density 
correlation function which in turn leads to a power-law dependence of the mean "mass" m(i?) within a domain of 
radius R around xq — (J: 

m{R)r^R°f, (9) 

with the exponent D/ in |(5J being the so-called fractal mass dimension P. [T^[r^ . Alternatively, using the saddle-point 
ensemble defined above the mean "mass" can be calculated as: 



m(i?) = (/ MxMO)\d''x). 

JR 



(10) 
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For a critical system at thermal equilibrium -D/ is related to the critical exponent 5, appearing in the effective action 
^ , as well as the dimension D of the embedding space through : 

The analysis at the level of one cluster is not sufficient for the description of the properties of the entire critical 
system. Therefore, after obtaining the set of saddle-point solutions of the generalized action JH)), it is necessary to use 
it in order to generate an ensemble of configurations valid for the global system. The systematic way to perform this 
extension is explained in detail in the next section. 



III. GENERATION OF MULTI-CLUSTER CONFIGURATIONS FOR THE CRITICAL SYSTEM 

In order to generate an ensemble of global configurations for the critical system within the saddle point approach 
discussed in the previous section one adopts the following picture: The entire system is composed of weakly correlated 
clusters of different size up to the correlation length |13J (which at the critical point is expected to be of the order of 
the lattice size). Thus, for given number and size of clusters, the partition function of the system can be decomposed 
in the product of the partition functions of each cluster. The sum over the microstates in the partition function of 
one cluster with given size is calculated through the summation over the corresponding saddle points. To saturate 
the summation over microstates in the total partition function one has to sum over the possible number of clusters 
as well as the corresponding possible cluster sizes. Therefore, the functional measure of the global partition function 
can be expressed as: 

P[^]^-r[.,vi ^Y^E^Iif ^^k^le-^I^-'^'l, (12) 
M {n,} ■ i=l 

where uflj^fii — Qt and nflj — ^ with i ^ j, in order to ensure that the M considered clusters of volume 
r^i are uncorrelated and cover the entire system of volume J^t- The functional integration J is defined as the 

sum over the saddle point solutions within the i-th cluster. According to the discussion in section II the dominant 
contribution to this sum comes from constant field configurations inside each cluster. Thus: 

/oo 
-oo 

In subsection (III A) we develop an algorithm for the integration in H12|) based on the saddle point approach in 1 — 13 
discussed so far. The algorithm is extended for higher dimensional systems in subsection (III B). Emphasis is given 
in the determination of observables carrying the fractal properties of the critical system. 



A. One Dimensional Case 



Our aim is to produce an ensemble of configurations, leading to a fractal mass dimension dictated by the power- 
law form of the density-density correlation of the critical system, in realistic computational times. According to the 
saddle point approach, a suitable basis to express the configurations of the critical system are the piecewise constant 
configurations, where the domains of constant value of the field ip{x) correspond to different clusters. In fact this 
coarse- graining procedure clearly neglects any internal structure of the critical clusters. However the entire set of 
these configurations, weighted with the corresponding Boltzmann factor calculated from the effective action, allows 
for the definition of suitable observables at the level of ensemble averaged quantities, reflecting the fractal geometry 
of the critical system. The main observable used in the present work is the mean "mass" in a domain of radius R, 
centered at xq, defined in ec ljlUI) . Therefore, the calculated value of the mass dimension Dj = '""j^f^-* will serve as a 
consistency check for the validity of our procedure. To construct the multi-cluster piecewise constant configurations 
we firstly perform a random partitioning of the iV-site lattice in a random integer number of elementary clusters M of 
different length. This is achieved using the Random Partition of an Integer n (RanPar) algorithm of [13 ■ Moreover, 
we use the Random Permutation of n Letters (RanPer) algorithm from the same reference, in order to permute and 
randomly select one specific partition. We point out here that this step is one of the most time consuming of the 
whole procedure, especially for iV > 4 x 10^. In the end, we come up with a random partitioning of the lattice into 
several clusters, and each cluster consists of a different number of lattice sites. 
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Continuing, we give a random constant value of the field (using a uniform distribution in the interval 
[~^maxT +^Tnax]) to the lattice sites within each cluster, and for this piecewise constant configuration we calculate 
the 1 — D eff'ective action ^ as: 



using periodic boundary conditions. Obviously, the derivative is zero inside each cluster and becomes relevant only 
in their edges. The choice of ifmax, of the lattice spacing a and of the coupling g, will be discussed later. 

As a next step we randomly alter the field value of the first cluster and we recalculate the effective action Fqi . If the 
new effective action Fqi is smaller than the initial Too we accept the new field value of the first cluster. On the other 
hand, if Fqi > Too we compare e''"'"'"'""^^ with a uniformly distributed random number z in the interval [0, 1), and if 
g(roo-roi) > 2; we also accept the change, otherwise we reject it. In this way the field configuration is weighted by a 
factor of e^'"['^l. We repeat this step changing the field value of the second cluster and calculating the new effective 
action ro2, and so on, until we have covered each one of the M clusters. In the end of this Metropolis algorithm 
we come up with a new field configuration and the corresponding action Fi — Tqm- We repeat this procedure at least 
kiter times in order to achieve equilibrium. The equilibration time kuer is defined as the minimum number of steps 
required in order to reach a stationary state, in the sense that variations of (p for Us (« 0(10'^)) successive algorithmic 
steps lead to a standard deviation of the mean value of cp and (p'^ less than 0.5%. The first configuration fulfilling this 
criterion is actually the first member of the ensemble of critical configurations and gets recorded. 

Repetition of the whole procedure will form the complete critical ensemble of the field configurations. Our algorithm 
allows for the use of a specific lattice partition more than once, in order to save computational time, since the 
partitioning procedure, using RanPar algorithm, as well as the equilibration process are time-consuming. Doing so, the 
necessary minimum iteration number luer, separating in time two successive statistically independent configurations, 
is much smaller than kner (usually luer ~ ~Lo!f )• However, the configuration space of the system is more efficiently 
covered using different lattice partitions, and therefore an optimization is possible leading to the use of the same 
partition for 10-50 times for an ensemble of > 10^ configurations. 

Finally, a comment concerning the choice of the various computational parameters has to be made. Firstly, the 
minimum iteration number luer can be calculated as usual in terms of the autocorrelation function G{m): 



where pn is the spatial average of the field in a configuration obtained by a single Metropolis step and time averaging 
for a large number of Metropolis steps (> 10"*) is performed. If m* is the characteristic decay time of G{m) then we 
choose liter ~ 5m*. 

Secondly, we have to find Pmax, which determines the optimal range of p values to be used in the simulation. 
Fortunately, there is a physical reason which constrains ipmax- Due to the potential term in the effective action |(5Jl, 
the statistical weight of configurations involving large if values are super-exponentially suppressed. Thus in thermo- 
dynamic equilibrium only configurations with field values within a narrow interval around zero will be statistically 
significant. 

For every choice of ipmax we produce a large number (~ 10**) of configurations and from this ensemble we calculate 
the average of the absolute value of the field (IfDi for every lattice site i. The quantity {\p\)i versus i corresponding 
to each Pmax is almost constant (with an expected noise that decreases increasing the ensemble population), with 
two peaks in the two lattice ends. While Pmax increases, a shift to larger constant {\(p\)i values is observed. However, 
above a threshold value for (pmax a saturation takes place. This behavior is expected since for small (fmax statistically 
significant field values contributing to the partition function are left out of the corresponding sum, while for Pmax^s 
larger than a specific value, the additional field values have vanishingly small contribution suppressed by the weight 
g-r[</5]^ So fmax can be strictly determined using a lower limit in the variation of the (IpDi vs i dependence. Knowing 
P'max reduces substantially the computational time as the requisite luer increases rapidly with increasing Pmax- 

Lastly, we have to fix the coupling g and the lattice spacing a. The former can be determined by considering the 
1 — D effective action as the projection of (O in one dimension, since in a genuine 1 — D system there is no critical 
point in the absence of long range interactions |l2l |. However, technically, an ensemble of configurations reproducing 
equivalently the density-density correlation of a fractal set with prescribed fractal mass dimension {Df < 1) can be 
constructed although having no direct physical interpretation. Concerning the lattice spacing a we have considered 
the thermodynamical and the continuum limit investigating the behavior of the term (^) for — > oo and a = const 
or a — > 0. In both cases we have observed a smooth behavior, attributed to the used basis of piecewise constant 
configurations, as well as to the the effect of the statistical weight e~^^'^^ of each configuration which favors smoother 
ones as TV — !■ oo. Therefore, the proposed algorithm is in fact independent of a. However, if one desires to evolve the 




(13) 



G(m) = {ipHPl + m) ~ {ipi){ipi+m), 



(14) 
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produced configurations in time, it is necessary to use the same lattice spacing in the numerical evolution, in order to 
achieve self-consistency 10]. Finally, we stress that in our presented results the calculated observables are expressed 
in units of a. 

With the demonstrated procedure we acquire an ensemble of field configurations generating a random fractal 
measure on the lattice as a statistical property after ensemble averaging. This fractal property lO is not directly 
reflected in the geometry of their average {(p{x)), but is produced only through the entire ensemble. The fractal mass 
dimension is determined by the power-law behavior of 



m{R) = {i \^{x)^mdx) 



(15) 



around = 0, averaged inside clusters of size R. The (J^ |(/7(a;)(^(0)| dx) versus R figure is drawn as follows: For 
a = of a specific configuration we find the size R of the cluster in which it belongs and we calculate the 
integral ^ ^\'^{x)ip{Q)\dx, thus acquiring one point in the {j ^\ip{x)'^{Q)\dx) vs R figure. For the same a;o = we 
repeat this procedure until we cover the whole ensemble, and the aforementioned figure is formed. Taking a different 
reference point xo obviously does not alter the results, since due to translation invariance in the averaged quantities 
m{xQ, R) ~ m[xQ + I, R), with I spanning the entire lattice. Eventually, our algorithm generates an ensemble of field 
configurations with fractal mass dimension 



Df 



(16) 



As an application we produce an ensemble of 10'' one dimensional field configurations on a = 2000 lattice, using 
(5 = 5 and g = 2, in which case the theoretical value of the fractal mass dimension according to (|16|l is 5/6. In fig. |21 
we depict the ensemble average {f{x)) having a noisy profile However, in the log-log plot of (J^ |(p(x)93(0)| dx) 
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FIG. 2; The (p-field on the 1 — D lattice averaged over the ensemble of the imtial configurations. 

versus R depicted in fig. 01 the slope, i.e the fractal mass dimension Df according to lfTC|l . is equal to 5/6 within an 
error of less than 0.3%. 

B. Higher Dimensional Case 



The higher dimensional generalization of our algorithm is straightforward, preserving the improved mean field 
approach using piecewise constant configurations. For D > 1 we use configurations consisting of D— dimensional 
boxes as basis and the ensemble production is reduced to the Cartesian product of the one dimensional case. Finally, 
the decisive test about the proper configuration production will be the calculation of the fractal mass dimension H10|l , 
as in the 1 — D case. In particular the 3 — D case is physically very interesting since the effective action (|2Jl describes 
the order parameter dynamics of the Ising universality class at the spontaneous magnetization transition point. Thus, 
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FIG. 3: (J^ \ip{x)Lp{0)\ dx) versus R for the ensemble of ip- field configurations. The slope, i.e the fractal mass dimension Df, 
is equal to 5/6 within an error of less than 0.3%. 



contrary to the one dimensional case, the investigated 3 — D system has a physical impact describing the standard 
model of a very common critical state. 

We firstly perform a random partitioning of the D— dimensional lattice in a random integer number of elementary 
box-shaped clusters of different volume, each one consisting of several lattice points. This is succeeded by applying 
the 1 — D partitioning algorithm D times and then taking the Cartesian product of the 1~ D partitions. In this case 
the linear size of the lattice is reduced, leading to a significant decrease of the partitioning algorithm computational 
time. For example, the time needed for the partitioning of a 20 x 20 x 20 lattice is two orders of magnitude smaller 
than that of a 2000— site linear lattice. 

Similarly to the \ — D case, we assign a random constant value of the field (in the interval [—ipmax, +^max]) to the 
lattice sites within each cluster, and for this piecewise constant configuration we calculate the Z?-dimensional effective 
action |(2J), using the corresponding generalized formula of eq. (I13f) . The sum extends to every lattice site using periodic 
boundary conditions, and the derivative is calculated using the straightforward D-dimensional generalization. The 
value of the coupling g in the 3 — D Ising effective action has been determined in the literature {g k, 2 4, 5j ) and 
the parameters ^Pmax-^ ^iter and luer are determined following the corresponding \ — D steps. The only complication 
enters in the specification of ipmax, which requires the construction of the plot {\(p\)i versus lattice site i, as we have 
mentioned in the previous subsection, becoming computationally much more demanding in this higher dimensional 
case. 

As a next step, we randomly change the field value of the first box-shaped cluster, we recalculate the effective 
action and using the same criteria as in the 1 — D case either we accept or reject the new field value of the first 
cluster. We perform these steps until we have covered the whole lattice, and we iterate this procedure kner times, in 
the end of which we record one field configuration. Repetition of the above algorithm will form the whole ensemble 
of the field configurations. Finally, our comment in the 1 — D treatment about the multiple use of a specific lattice 
partition, holds in the present case, too. However, due to the significantly smaller partitioning computational time, 
such a treatment is not necessary. 

This is the generalized algorithm of producing an ensemble of I?-dimensional field configurations generating a 
random fractal measure on the lattice. The ensemble possesses the property of eqs. IjlOII and where now the 
fractal mass dimension is determined by the power-law behavior of m(R) around xq = 0, averaged inside clusters of 
volume V. Similarly to the 1 — D case, to acquire the (J^ \(p{x)ip{0) \ d^x) versus R figure we chose xq = of a specific 
configuration, we find R of the cluster to which it belongs, approximated as cx \/V ^ and we calculate the integral 
/ \lp{x)lp{u)\ d^x, thus obtaining one point in the {Jj^ \(p{x)(p{0)\ d^x) vs R figure. Repetition of this procedure for 
the whole configuration ensemble provides all the points in Fig. ^ 

As an application we produce an ensemble of 10* 3-dimensional field configurations on a 20 x 20 x 20 lattice, using 
5 — b and g = 2 As already discussed, this choice has a physical correspondence, since for dimensionality 

D — ?>, isothermal critical exponent (5 = 5 and coupling g — 2, the free energy describes the effective action of the 
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3 — -D Ising model at its critical point. In this case, the theoretical value of the fractal mass dimension according to 
(HU is 5/2. 

In Fig. 0] we observe that in the log-log plot of (J^ |(/3(x)(^(0)| cPx) vs R, the slope, i.e the fractal mass dimension 
Df according to Q, is equal to 5/2, within an error of less than 1%. In order to test that the proposed algorithm is 




10 

R (in units of a) 



FIG. 4: { \(p{x)ip{0) I d^x) versus R for the ensemble of3 — D ip-field configurations. The slope, i. e the fractal mass dimension 
Df, is equal to 5/2 within an error of less than 1%. 

valid for any lattice site, as dictated by the statistics in the considered system, we have calculated Df using different 
reference points xq (sources in on the 3 — D lattice. The corresponding distribution p{Df) is shown in Fig. |S1 
It is clearly seen that Df is almost constant with a deviation of at most 4% for the given lattice size N^, and as we 



1500- 




FIG. 5: p{Df) calculated using 8 x 10^ different reference points xq on the lattice. 

have tested, increasing N this deviation decreases algebraically. 

Finally, a last test is performed in order to check the ability of the obtained ensemble of configurations to reproduce 
the statistical properties of the critical system. Besides the underlying fractal geometry of the critical clusters, the 
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two-point correlation function: 



C{x,y)^Mx)^{y)\)-Mx)\)My)\) 



(17) 



possesses an analytically known power-law form at the critical point. Using the constructed ensemble of configurations 
we have calculated numerically the correlation function H17|l and the result is shown in Fig. |H| The theoretical 
expectation C{x,y) ^ \x — y\~^~^ (?? « 0.04 for the 3 — D Ising universality class), corrected by a small exponential 
factor incorporating finite size effects, is shown with the dashed line. It is clearly seen that the calculated C{x, y) is 
in very good agreement with the analytical formula supporting further the equivalence of the obtained ensemble of 
configurations with the critical state of the considered system. 
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FIG. 6; C{x, y) as a function of \x — y\ calculated using 10* critical configurations for the 3 - 
(Ising). The dashed line presents the theoretically expected result. 



D one-component real scalar field 



IV. CONCLUDING REMARKS 

Using the saddle-point approach introduced in we have been able to develop an algorithm simulating the 
critical state of a macroscopic system at its transition point. The method followed here is in close analogy with the 
improved mean field theory (2H and leads to a successful and computationally efhcient description of the geometrical 
characteristics of the critical clusters, defining suitable measures to quantify this property. A particularly appealing 
issue of the proposed method is that it provides a link between the effective action at the critical point and the 
fractal geometry of the formed clusters, overcoming the huge numerical effort needed for the detailed description 
of fractal sets. Thus, it is the first time, at least to our knowledge, that geometrical characteristics as well as 
statistical properties of the critical system, are incorporated in an ensemble of configurations suitable for the study 
of any desired observable at the critical point. The proposed algorithm may be of special interest for the treatment 
of systems where the formed critical state is not observable but acts as an intermediate state for the subsequent 
evolution of the system. Such a scenario is likely to hold in the collision of two heavy nuclei at high energies 
when the formed fireball freezes out near the theoretically predicted QCD critical point. 
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